# import numpy as np
# import netCDF4 as nc
# import matplotlib.pyplot as plt
# import cartopy.crs as ccrs
# import cartopy.feature as cfeature
#
# # 加载数据
# nf = nc.Dataset(r'/mnt/datastore/liudddata/result/20200101test/2020010104_predicted_2d.nc', 'r')
# cbh_data = np.ma.getdata(nf.variables['predicted'][:])  # 转换为普通数组
# lat = nf.variables['lat'][:]
# lon = nf.variables['lon'][:]
# nf.close()
#
# # 处理无效值（假设填充值为-999或其他特殊值）
# cbh_data = np.ma.masked_invalid(cbh_data)  # 自动屏蔽NaN和inf
# cbh_data = np.clip(cbh_data, 0, None)      # 去掉负值
#
# # 创建地图画布
# fig = plt.figure(figsize=(12, 8))
# ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
#
# # 添加地理特征
# ax.add_feature(cfeature.COASTLINE)
# ax.add_feature(cfeature.BORDERS, linestyle=':')
# ax.add_feature(cfeature.LAND, facecolor='lightgray')
#
# # 使用contourf绘制云高分布
# levels = np.linspace(0, cbh_data.max(), 20)  # 自动设置20个色阶
# contour = ax.contourf(lon, lat, cbh_data,
#                       levels=levels,
#                       cmap='Blues',
#                       extend='both',        # 扩展颜色条箭头
#                       transform=ccrs.PlateCarree())
#
# # 添加颜色条
# cbar = plt.colorbar(contour, orientation='horizontal', pad=0.05, aspect=30)
# cbar.set_label('Cloud Base Height (meters)')
#
# # 设置标题和范围
# ax.set_title('Cloud Base Height Distribution (2020-01-01 04UTC)')
# ax.set_global()  # 显示全球范围
#
# # 保存或显示图像
# plt.savefig('cbh_20200104.png', dpi=500, bbox_inches='tight')
# plt.show()

"""修改投影方式"""

import numpy as np
import netCDF4 as nc
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

# 加载数据
nf = nc.Dataset(r'/mnt/datastore/liudddata/result/20200104new/2020010104_predicted_2d.nc', 'r')
cbh_data = np.ma.getdata(nf.variables['predicted'][:])  # 转换为普通数组
lat = nf.variables['lat'][:]
lon = nf.variables['lon'][:]
nf.close()

# 处理无效值（假设填充值为-999或其他特殊值）
cbh_data = np.ma.masked_invalid(cbh_data)  # 自动屏蔽NaN和inf
cbh_data = np.clip(cbh_data, 0, None)      # 去掉负值

# 创建地图画布，使用正交投影
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.Orthographic(central_latitude=0, central_longitude=104.7))

# 添加地理特征
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle=':')
# ax.add_feature(cfeature.LAND, facecolor='lightgray') #空白区域为灰色
ax.add_feature(cfeature.LAND, facecolor='white')

# 添加经纬度线
parallels = range(-90, 91, 30)  # 每30度一个纬线
meridians = range(-180, 181, 30)  # 每30度一个经线
gl = ax.gridlines(draw_labels=True, color='gray', linestyle='--', xlocs=meridians, ylocs=parallels)
gl.top_labels = False
gl.right_labels = False
gl.xlabel_style = {'size': 10}
gl.ylabel_style = {'size': 10}

# 使用contourf绘制云高分布
levels = np.linspace(0, cbh_data.max(), 10)  # 自动设置20个色阶
contour = ax.contourf(lon, lat, cbh_data,
                      levels=levels,
                      cmap='Blues',
                      extend='both',        # 扩展颜色条箭头
                      transform=ccrs.PlateCarree())

# 添加颜色条
cbar = plt.colorbar(contour, orientation='vertical', pad=0.01, aspect=20)
cbar.set_label('Cloud Base Height (meters)')

# 设置标题
ax.set_title('Cloud Base Height Distribution (2020-01-01 04UTC)')

# 保存或显示图像
plt.savefig('cbh_20200104.png', dpi=300, bbox_inches='tight')
plt.show()


'''平均云底高度'''
# import numpy as np
# import netCDF4 as nc
# import matplotlib.pyplot as plt
# import cartopy.crs as ccrs
# import cartopy.feature as cfeature
# import os
# import glob
#
# # 设置数据路径
# data_dir = "/mnt/datastore/liudddata/result/20200101test"
# output_dir = "./results"
# os.makedirs(output_dir, exist_ok=True)
#
# # 获取所有NC文件列表
# nc_files = sorted(glob.glob(os.path.join(data_dir, "*_predicted_2d.nc")))
# print(f"找到 {len(nc_files)} 个数据文件")
#
# # 初始化数据存储列表
# all_cbh = []
#
# # 读取所有文件数据
# for i, file in enumerate(nc_files):
#     try:
#         with nc.Dataset(file, 'r') as nf:
#             # 读取数据并转换为普通数组
#             cbh = np.ma.getdata(nf['predicted'][:])
#             lat = nf['lat'][:]
#             lon = nf['lon'][:]
#
#             # 处理无效值（假设使用-999作为填充值）
#             cbh = np.ma.masked_where(cbh < 0, cbh)  # 过滤负值
#             cbh = np.ma.filled(cbh, np.nan)  # 将masked值转为nan
#
#             all_cbh.append(cbh)
#
#         print(f"已加载 {i + 1}/{len(nc_files)} 文件: {os.path.basename(file)}")
#     except Exception as e:
#         print(f"错误读取文件 {file}: {str(e)}")
#         continue
#
# # 转换为三维数组 [时间, 纬度, 经度]
# all_cbh = np.array(all_cbh)
#
# # 计算时间平均（忽略nan）
# mean_cbh = np.nanmean(all_cbh, axis=0)
#
# # 创建绘图画布
# fig = plt.figure(figsize=(12, 8))
# ax = fig.add_subplot(111, projection=ccrs.PlateCarree())
#
# # 添加地理特征
# ax.add_feature(cfeature.COASTLINE, linewidth=0.8)
# ax.add_feature(cfeature.BORDERS, linestyle=':', linewidth=0.5)
# # ax.add_feature(cfeature.LAND, facecolor='lightgray')
# ax.add_feature(cfeature.LAND, facecolor='white')
# # ax.add_feature(cfeature.OCEAN, facecolor='lightcyan')
# ax.add_feature(cfeature.OCEAN, facecolor='white')
#
# # 设置色阶参数
# vmin = np.nanpercentile(mean_cbh, 5)  # 取5%分位数作为最小值
# vmax = np.nanpercentile(mean_cbh, 95)  # 取95%分位数作为最大值
# levels = np.linspace(vmin, vmax, 15)  # 生成15个等间距色阶
#
# # # 绘制等值线填充图
# # cf = ax.contourf(lon, lat, mean_cbh,
# #                  levels=levels,
# #                  cmap='YlGnBu',
# #                  extend='both',
# #                  transform=ccrs.PlateCarree())
# # 绘制等值线填充图
# cf = ax.contourf(lon, lat, mean_cbh,
#                  levels=levels,
#                  # cmap='YlGnBu',
#                 cmap='Blues',
#                  extend='both',
#                  transform=ccrs.PlateCarree())
#
# # 添加颜色条
# cbar = plt.colorbar(cf, orientation='horizontal',
#                     pad=0.05, aspect=30,
#                     format='%.0f')
# cbar.set_label('Cloud Base Height (m)', fontsize=12)
#
# # 添加标题
# plt.title(f'Average Cloud Base Height ({len(nc_files)} files average)',
#           fontsize=14, pad=20)
#
# # # 设置地图范围（示例设置中国区域）
# # ax.set_extent([70, 140, 15, 55], crs=ccrs.PlateCarree())
# # 设置标题和范围
# ax.set_title('Cloud Base Height Distribution (2020-01-01 00UTC)')
# ax.set_global()  # 显示全球范围
#
# # 保存图片
# output_path = os.path.join(output_dir, 'average_cbh.png')
# plt.savefig(output_path, dpi=500, bbox_inches='tight')
# print(f"图片已保存至 {output_path}")
#
# # 显示图像
# plt.show()